import numpy as np
import matplotlib.pyplot as plt
import math

def calculate_e_with_theorems():
    """
    使用夹逼定理和单调有界定理计算自然常数e
    """
    # 设置中文字体
    plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
    plt.rcParams['axes.unicode_minus'] = False
    
    # 创建图形
    fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
    
    # 参数设置
    n_max = 100
    n_values = np.arange(1, n_max + 1)
    
    # 定义两个关键数列
    # a_n = (1 + 1/n)^n 单调递增
    # b_n = (1 + 1/n)^(n+1) 单调递减
    a_n = (1 + 1/n_values)**n_values
    b_n = (1 + 1/n_values)**(n_values + 1)
    
    # 真实e值
    e_true = math.e
    
    print("=" * 60)
    print("使用夹逼定理和单调有界定理计算自然常数e")
    print("=" * 60)
    
    # 验证夹逼定理条件
    print("\n1. 夹逼定理验证：")
    print("   - 夹条件：a_n ≤ e ≤ b_n 对所有n成立")
    print("   - 逼条件：lim(n→∞) a_n = lim(n→∞) b_n = e")
    print("   因此，根据夹逼定理，e = lim(n→∞) a_n = lim(n→∞) b_n")
    
    # 验证单调有界定理条件
    print("\n2. 单调有界定理验证：")
    print("   - a_n 单调递增：a_{n+1} ≥ a_n")
    print("   - b_n 单调递减：b_{n+1} ≤ b_n") 
    print("   - a_n 有上界：a_n ≤ b_n ≤ 4")
    print("   - b_n 有下界：b_n ≥ a_n ≥ 2")
    print("   根据单调有界定理，两个数列均收敛")
    
    # 在第一个子图上绘制夹逼定理
    ax1.plot(n_values, a_n, 'b-', label=r'$a_n = (1 + \frac{1}{n})^n$ (递增)', linewidth=2)
    ax1.plot(n_values, b_n, 'r-', label=r'$b_n = (1 + \frac{1}{n})^{n+1}$ (递减)', linewidth=2)
    ax1.axhline(y=e_true, color='green', linestyle='--', 
                label=f'真实值 e ≈ {e_true:.10f}', linewidth=2)
    
    # 标记夹逼区间
    for n in [1, 5, 10, 20]:
        ax1.plot([n, n], [a_n[n-1], b_n[n-1]], 'k-', alpha=0.3)
    
    ax1.set_title('夹逼定理：a_n ≤ e ≤ b_n')
    ax1.set_xlabel('n')
    ax1.set_ylabel('函数值')
    ax1.grid(True, linestyle='--', alpha=0.7)
    ax1.legend()
    ax1.set_ylim(2.5, 3.0)
    
    # 在第二个子图上展示误差收敛
    a_error = np.abs(a_n - e_true)
    b_error = np.abs(b_n - e_true)
    interval_width = b_n - a_n
    
    ax2.semilogy(n_values, a_error, 'b-', label=r'$|a_n - e|$', linewidth=2)
    ax2.semilogy(n_values, b_error, 'r-', label=r'$|b_n - e|$', linewidth=2)
    ax2.semilogy(n_values, interval_width, 'g-', label='夹逼区间宽度', linewidth=2)
    
    ax2.set_title('误差分析和收敛速度')
    ax2.set_xlabel('n')
    ax2.set_ylabel('误差（对数坐标）')
    ax2.grid(True, linestyle='--', alpha=0.7)
    ax2.legend()
    
    plt.tight_layout()
    plt.show()
    
    # 数值计算结果
    print("\n3. 数值计算结果：")
    print("-" * 60)
    print(f"真实值 e = {e_true:.15f}")
    print("\nn\t a_n\t\t\t b_n\t\t\t 区间宽度\t a_n误差")
    print("-" * 90)
    
    test_n = [1, 5, 10, 50, 100, 1000]
    for n in test_n:
        a_val = (1 + 1/n)**n
        b_val = (1 + 1/n)**(n+1)
        width = b_val - a_val
        a_err = abs(a_val - e_true)
        
        print(f"{n}\t {a_val:.10f}\t {b_val:.10f}\t {width:.2e}\t {a_err:.2e}")
    
    # 验证单调性
    print("\n4. 单调性验证：")
    print("-" * 60)
    
    # 验证a_n单调递增
    a_monotonic = all(a_n[i] <= a_n[i+1] for i in range(len(a_n)-1))
    b_monotonic = all(b_n[i] >= b_n[i+1] for i in range(len(b_n)-1))
    
    print(f"a_n 单调递增: {a_monotonic}")
    print(f"b_n 单调递减: {b_monotonic}")
    
    # 验证夹逼条件
    squeeze_condition = all(a_n[i] <= e_true <= b_n[i] for i in range(len(a_n)))
    print(f"夹逼条件满足 (a_n ≤ e ≤ b_n): {squeeze_condition}")
    
    # 计算极限近似值
    print("\n5. 极限计算结果：")
    print("-" * 60)
    
    # 使用较大的n计算极限近似值
    large_n = 1000000
    e_approx = (1 + 1/large_n)**large_n
    error = abs(e_approx - e_true)
    
    print(f"使用 n = {large_n} 计算：")
    print(f"近似值: {e_approx:.15f}")
    print(f"真实值: {e_true:.15f}")
    print(f"绝对误差: {error:.2e}")
    print(f"相对误差: {error/e_true:.2%}")

def theoretical_proof():
    """
    定理的数学证明概要
    """
    print("\n" + "=" * 60)
    print("数学证明概要")
    print("=" * 60)
    
    print("""
1. 单调有界定理证明：
   - 使用二项式定理展开 a_n = (1 + 1/n)^n
   - 证明每一项都是非负且递增的
   - 证明数列有上界（如 3）
   - 因此极限存在

2. 夹逼定理证明：
   - 构造不等式：a_n ≤ e ≤ b_n
   - 证明 lim(n→∞) (b_n - a_n) = 0
   - 因此 lim(n→∞) a_n = lim(n→∞) b_n = e

3. 极限计算：
   e = lim(n→∞) (1 + 1/n)^n
     = 1 + 1/1! + 1/2! + 1/3! + ... (泰勒级数展开)
    """)

if __name__ == "__main__":
    calculate_e_with_theorems()
    theoretical_proof()
